#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#include <sys/times.h>
#include "mytime.h"

#define N_DEF 200
#define TOL_DEF 0.0000001
#define MAXIT_DEF 500

/****************/
int *despl;
int *count;
/****************/

double jacobi_sec(double **a, int n, double *b, double *x, double tol, int maxit) {
  int i, j, it = 0;
  double sum, resi = 10.0, *new_x;

  new_x = (double *) malloc(n*sizeof(double));

  for (i = 0; i < n; i++) x[i] = b[i];
  /* while ( (it < maxit) && (resi > tol) ) { */
  while (it < maxit) {
    resi = 0.0;
    for (i = 0; i < n; i++) {
      sum = 0.0;
      for (j = 0; j < n; j++) 
        sum += a[i][j]*x[j];
      resi += fabs(sum - b[i]);
      sum += -a[i][i]*x[i];
      new_x[i] = (b[i] - sum) / a[i][i];
    }
    for (i = 0; i < n; i++) x[i] = new_x[i];
    it++;
  } 

  free(new_x);
  return resi;  
}

int main(int argc, char *argv[]) {
  double tol;
  double **a, *b, *x;
  int i, j, n, diag, maxit;

  CLOCK_TYPE chrono;
  double seq_time;
  double seq_resi;

  n = N_DEF;
  maxit = MAXIT_DEF;
  tol = TOL_DEF;

  switch (argc) {
    case 4: tol = atof(argv[3]);
    case 3: maxit = atoi(argv[2]);
    case 2: n = atoi(argv[1]);
  }

  ReservaMatriz(&a, n, n);
  ReservaVector(&b, n);
  ReservaVector(&x, n);

    GeneraVector(b, n);

  /* matriz aleatoria */
     GeneraMatriz(a, n, n);
     /* Para garantizar la convergencia del metodo de Jacobi */
     /* Convierte en diagonal dominante la matriz distribuida */
     diag = 0;
     for (i = 0; i < n; i++) {
        for (j = 0; j < n; j++) 
           a[i][diag] += a[i][j];
        diag++;
     }
 
  /****************************************************************************
   algoritmo secuencial */
  CLOCK_Start(chrono);
  seq_resi = jacobi_sec(a, n, b, x, tol, maxit);
  CLOCK_End(chrono, seq_time);
 /* CompruebaSolucion (a, x, b, n); */

  printf("%d:%g:%g\n", 0, seq_time, seq_resi); 

  
  /**************************************************************************
   Libeacion de memoria */
  LiberaMatriz(a);
  LiberaVector(b);
  LiberaVector(x);

  return 0;

}
